function IdSet=OneSide(u0grid, u1grid, u2grid, u3grid, ...
                       sgreduced, PYX, ...
                       id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
                       P_col,...
                       Tcoord_zeros,...
                       Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
                       cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros, equalorder, n_vertices,infty_pos,...
                       cr_symmetry, fill_symmetry, RHS_symmetry, S1, S2, S3)

exists=zeros(sgreduced,1);
param_MOSEK.MSK_IPAR_LOG = 0; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Outside the loop across parameter values %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% RHS 
RHS_obsequiv=(-1+PYX).';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
%%%%%%%%%%%%%%%% Inside the loop across parameter values %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prob=[];

parfor g=1:sgreduced
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:); 
idU3=id_Ugridreduced{3}(g,:); %for each of the U-sets
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
%s3=s_id_gridreduced(g,3);
sU=sU_id_gridreduced(g); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Equalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% Coordinate columns
cc_obsequiv=[changecoord_3(s1,s2, idU1(P_col{1}(:,1)), idU2(P_col{1}(:,2)), idU3(P_col{1}(:,3))) ... %Y=1
             changecoord_3(s1,s2, idU1(P_col{2}(:,1)), idU2(P_col{2}(:,2)), idU3(P_col{2}(:,3))) ... %Y=2
             changecoord_3(s1,s2, idU1(P_col{3}(:,1)), idU2(P_col{3}(:,2)), idU3(P_col{3}(:,3))) ... %Y=3
             changecoord_3(s1,s2, idU1(P_col{4}(:,1)), idU2(P_col{4}(:,2)), idU3(P_col{4}(:,3)))];   %Y=0 

%% Integrates to one
cc_one=changecoord_3(s1,s2, idU1(infty_pos), idU2(infty_pos), idU3(infty_pos));

%% Zeros
cc_zeros=changecoord_3(s1,s2, idU1(Tcoord_zeros(:,1)), idU2(Tcoord_zeros(:,2)), idU3(Tcoord_zeros(:,3))); %1xn_zeros

%% Symmetry
cc_symmetry=[reshape([changecoord_3(s1, s2, idU1(S1(:,1)), idU2(S1(:,2)), idU3(S1(:,3))).' ... 
                      changecoord_3(s1, s2, idU1(S2(:,1)), idU2(S2(:,2)), idU3(S2(:,3))).'].', 1, (n_U_sets+1)*n_U_sets*2) ...
             ...         
             changecoord_3(s1,s2, idU1(S3(:,1)), idU2(S3(:,2)), idU3(S3(:,3)))]; 
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Inequalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Increasingness 
%Repeat the way we constructed Tcoord but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets

D=[T(indices_pairs_extended(:,1),:) T(indices_pairs_extended(:,2),:) ...
   Tcoord(indices_pairs_extended(:,1),:) Tcoord(indices_pairs_extended(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates


%Save in D1 the rows of D where the first n_U_sets-tuple is <= than the second
%Only the columns reporting column-coordinates are relevant
D1=D(sum(D(:,1:n_U_sets)<=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd1x(2*n_U_sets)
sd1=size(D1,1);

%Save in D2 the rows of D where the first n_U_sets-tuple is >= than the second
%Only the columns 5:end (reporting column-coordinates) are relevant
D2=D(sum(D(:,1:n_U_sets)>=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd2x(2*n_U_sets)
sd2=size(D2,1);

%Flip D2
D2=[D2(:,n_U_sets+1:2*n_U_sets) D2(:,1:n_U_sets)]; %first n_U_sets-tuple is <= than the second

%Combine D1 and D2
D=[D1;D2];
sd=sd1+sd2; 

cr_increas=kron((1:1:sd), ones(1,n_vertices)); %1x[sd*n_vertices]
%n_vertices comes from the number of all possible vertices generated by each pair of n_U_sets-tuples

cc_increas=zeros(1,sd*n_vertices);
for j=1:sd
    D_temp=D(j,:);
    vertices=D_temp(C); %n_verticesxn_U_sets
    cc_increas(n_vertices*(j-1)+1:n_vertices*j)=changecoord_3(s1,s2, idU1(vertices(:,1)), idU2(vertices(:,2)), idU3(vertices(:,3))); %1xn_vertices
end

fill_increas=repmat(-C_sign, 1, sd);  %1x[sd*n_vertices]
%remember: signs are flipped because inequalities are <= 

RHS_increas=zeros(sd, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Compose %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aeq=sparse([cr_obsequiv cr_zeros cr_one cr_symmetry], [cc_obsequiv cc_zeros cc_one cc_symmetry], [fill_obsequiv fill_zeros fill_one fill_symmetry], numeqconstraints, sU);   %[#equalityconstraints]x[sU]
Aineq=sparse(cr_increas, cc_increas, fill_increas, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]
beq=[RHS_obsequiv; RHS_zeros;RHS_one; RHS_symmetry];

% MOSEK
problocal=prob;
problocal.blx=zeros(sU,1); %lower bound unknowns
problocal.ulx=ones(sU,1); %upper bound unknowns
problocal.c=zeros(sU,1); %objective function
problocal.a=[Aeq; Aineq];
problocal.blc=[beq; -Inf(size(Aineq,1),1)];
problocal.buc=[beq; RHS_increas];
[~,result]     = mosekopt('minimize echo(0)',problocal, param_MOSEK);
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   exists(g)=1; 
end
end

%Enlarge results to the original grid of parameter vectors
exists=exists(equalorder);
IdSet = [u0grid(logical(exists),:)  u1grid(logical(exists),:)  u2grid(logical(exists),:) u3grid(logical(exists),:)];

end

